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Abstract 

We apply a multiresolution community detection algorithm to perform unsupervised segmentation 
of complex intracellular signals derived using fluorescent dyes. In our earlier work [f], when applying 
our method to benchmarks, our algorithm was shown to be one of the best and to be especially suited 
to the detection of camouflage images. In the current manuscript, we have explored this algorithm 
in a more complex scenario. The current image processing problem is framed as identifying clusters 
with respective average fluorescent lifetimes (FLTs) against a background or "solvent" in fluorescence 
lifetime imaging microscopy (FLIM) images derived using NIR fluorescent dyes. We have identified 
significant multiresolution structures using replica correlations in these images, where such correlations 
are manifested by information theoretic overlaps of the independent solutions ("replicas") attained using 
the proposed algorithm from different starting points. Our method is more efficient than a well-known 
image segmentation method based on mixture of Gaussian distributions. It offers more than 1.25 times 
diversity based on Shannon index than the latter method, in selecting clusters with distinct average FLTs 
in NIR FLIM images. 
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1 Introduction 



Image segmentation plays a crucial role in many medical imaging applications by enhancing the detection of 
anatomical structures of interest. Examples of medical image segmentation methods are listed in [2], which 
include graph partitioning methods and normalized cuts [3,4], and mixture of Gaussian distributions (MGD) 
method [5], to name a few. 

In this work, we propose a multiresolution "community detection" (CD) approach based on graph parti- 
tioning theory [3] to segment complex intracellular signals derived using fluorescent dyes. CD [6-11] seeks to 
divide groups of nodes with dense connections internally and with sparser connections between the groups. 
Moreover, it partitions a large physically interacting system into optimally decoupled communities. To 
demonstrate our approach, we have used fluorescence lifetime imaging microscopy (FLIM) data captured 
using near-infrared (NIR) fluorescent dyes [12], where the underlying signal describing intra-ccllular distri- 
bution is complex in nature. FLIM, a promising technique for imaging molecular process, uses time-resolved 
measurements of fluorescence from cells and thin tissue sections to generate images of the characteristic 
fluorescence lifetimes (FLTs) within a pixel or voxel. The FLT is the average time a molecule resides in the 
excited state before returning to the ground state through fluorescence emission [13]. 

Segmentation of the FLIM data is challenging due to the high amount of spatial and temporal noise 
attached to such data. Such a problem is severe for images captured using organic NIR dyes because of 
their low photostability, resulting in low-signal-to-noise ratio (SNR) and low-resolution images. Hence, there 
exists a niche to efficiently segment such data to mine the spatial structures hidden in them. In this study, 
we deliberately used low SNR data to determine the segmentation quality of our CD method. Although we 
focus on the NIR FLIM data in this report, the CD method is generally applicable to diverse imaging data. 

To perform the segmentation, our multiresolution CD first investigates the optimal structure at different 
resolutions in the input image data. It then analyzes each of these resolutions to obtain the respective 
number of estimated communities and the corresponding partition strengths. Partition strength is evaluated 
via information correlations between independent solutions ("replicas") starting from different initial states. 
This analysis determines the significant structures at which the replicas are strongly correlated. The outcome 
of the proposed method for the case of FLIM images is a segmented image containing distinct average FLTs 
in each of its segments. It is noteworthy to mention that our CD method is "unsupervised" in nature. 
Namely, it does not employ any ground-truth as a prior knowledge to train the algorithm. In particular, 
the algorithm does not assume any prior knowledge of specific spatial patterns corresponding to any specific 
community hidden in the input data. 

We have compared the performance of our proposed CD method with another image segmentation method 
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based on the mixture of Gaussian distributions (MGD) [5] for segmenting the NIR FLIM images. Our method 
proves to be more than 1.25 times diverse based on Shannon index than the latter method, in finding spatial 
regions with distinct average FLTs in the NIR FLIM data. A detailed biological validation of the resulting 
segments and their biological roles is beyond the scope of this work, since such validation varies from one 
biological study to another. We nevertheless argue that our proposed method is general for segmenting 
fluorescence microscopy images, and it can be used to generate hypotheses for future biological validation. 

The study is presented as follows. In Section 2, we discuss briefly on our NIR FLIM system and the 
structure of the data acquired from such system. Section 3 presents our proposed multiresolution CD 
approach for performing segmentation of images based on their spatial information. Section 4 extends our 
approach for images described by both spatial and temporal information. In both of these sections, we 
also develop image segmentation methods based on MGD. In Section 5, we present the performance of our 
proposed method using NIR FLIM image of an ex vivo liver tissue sample of mice treated with a NIR 
fluorescent dye. Section 5 describes a comparison of this performance with that attained using the MGD 
method. We conclude in Section 6. 

2 Near-Infared Fluorescence Lifetime Imaging Microscopy 

Before we present our segmentation method in detail, we briefly describe an overview of our NIR FLIM 
system and the structure of the acquired data [12]. Our NIR FLIM system consists of a fiber-coupled laser 
diode (BDL-785-SMC, Becker-Hickl, Gremany), a confocal laser scanning fluorescence microscope (FV1000, 
Olympus, Center Valley, PA), a thermoelectrically cooled, red-enhanced photomultiplier tube (PMT) (PMC- 
100-20, Becker-Hickl, Germany), and a time-correlated single photon counting (TCSPC) card (SPC-730, 
Becker-Hickl, Germany). The laser diode operates at TEM o mode [14] and provides 785 nm excitation 
light. Its pulsed wave duration is nominally 60-80 ps with frequency of 50 MHz. The PMT offers a minimum 
photon-count rate of 5 MHz. The TCSPC card has a transit-time spread of 180 ps and a dead time of 125 
ns. The single photon counting operation employs equally spaced 256 time gates of duration 16.6 ns, with 
an initial delay of 1.4 ns for eliminating photons from the excitation pulse in the measurement. 

To capture the FLIM images, the laser light was collimated, passed through the confocal system, and 
focused onto the sample using a 20X, 0.95-NA objective. Single photon fluorescence from the sample was 
collected through the same objective and directed by a dichroic mirror toward the confocal pinhole and 
detected by the PMT. Residual excitation light was removed using a bandstop filter with cutoff wavelengths 
of 765 nm and 805 nm before the PMT. Data acquisition was performed using the TCSPC card, triggered 
via a synchronization signal generated by the laser driver for each laser pulse. Images were acquired by 
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unidirectional scanning with the excitation beam using a galvanometric mirror pair as embedded in the 
confocal microscope system. 

To generate the input data for the segmentation methods used in this study, collected time traces captured 
from the NIR FLIM system were first analyzed in the SPCImage software (Becker- Hickl, Germany). It 
replaces time traces per pixel with the cumulative time traces computed using the target pixel and its 
neighboring eight pixels to increase the SNR. Maximum SNR here is typically defined by the square root 
of the number of photons acquired in the peak channel of the 256 time gates used to perform TCSPC 
acquisition [15]. The resulting data of size 128x128x256 describe spatio-temporal FLT information of the 
target field-of-view (FOV) of the imaging sample of interest. The first two dimensions of size 128x128 
describes the spatial pixel locations in this FOV. The third dimension of length 256 samples for any of 
these pixels describes a temporal convolution between the instrument response function (IRF) and the 
fluorescent decay trace of that particular pixel location. To generate a spatial two-dimensional (2D) dataset 
for performing the segmentation, the above-mentioned three-dimensional (3D) data were analyzed further 
using the SPCImage software. In each pixel location of the imaging sample, the earliest time points allow 
fitting of the IRF, while a single-exponential curve was sufficient to fit the subsequent falloff with time along 
the third dimension of length 256. A Xr fitness test determined the validity of the fit, providing x 2 r values 
< 1.5 for all pixels. The resulting data of size 128x128 describe a spatial FLT information of the target 
FOV of the imaging sample of interest. This report presents methods for performing spatial segmentation 
of both 2D and 3D FLIM data described herein. 

To conduct specific examples for analyzing performances of the segmentation methods explored in this 
work, we imaged ex vivo liver tissue sample of a mice treated with an NIR fluorescent dye. The NIR FLIM 
data was captured using the system and methodology as described above. All animal studies were performed 
in compliance with the Washington University School of Medicine Animal Studies Committee requirements 
for the humane care and use of laboratory animals in research. 

3 Segmentation of Images Described by 2D Spatial Information 

This section first discusses our proposed CD method for spatially segmenting images described by 2D spatial 
information. We also briefly review the MGD method for segmenting such images. 

3.1 Community Detection Method 

We employ the multiresolution CD algorithm to investigate the optimal structure at different resolutions 
in the input images. We use the number of estimated communities and the information based measures to 
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determine the significant structures at which the "replicas" are strongly correlated, ("replicas" are defined 
as independent solutions of the multiresolution CD algorithm attained from different starting points.) We 
determine different levels of detail and resolutions by setting the resolution parameters. 
We minimize a Potts model Hamiltonian for solving the CD problem [11], 



K 

SfatTi). (1) 

Here the number of communities (segments) K can be specified from the input or left arbitrary, allowing 
the algorithm to decide the number of segments K using the lowest energy solutions. The weight Wij denotes 
the absolute FLT difference between a pixel pair formed by the i th and j th € {1,2,..., TV}) pixels in 

the input image with N pixels. The Heavyside functions 0(-) "turns on" or "off" the edge designation. 



W) &(W - Wij) + jQ(W i:j - W) 



e(Wij -w) = { 



if > w, 

otherwise. 



(2) 



The Kronecker delta S(-) is given by Eq. (3), 



if Gi = <jj, 

\ (3) 
0, otherwise. 

In this Hamiltonian, by virtue of the 5(<Ti, Oj) term, each spin a interacts only with other spins in its own 
segment. The spin a, (Vcr^ G {1,2,..., K}) defines the segment identity for the i th (i & {1,2, ... , N}) pixel, 
and the algorithm optimizes it in minimizing the energy defined by Eq. (1). As such, the resulting model is 
local — a feature that enables high accuracy along with rapid convergence [10]. As before, minimizing the 
Hamiltonian of Eq. (1) corresponds to identifying strongly connected segments of pixels. 

For the 2D NIR FLIM images, we define the weight of any edge formed by two pixels as the absolute 
difference between their FLTs. We then apply the multiresolution CD algorithm to segment the network 
formed using the image pixels as nodes. We then analyze the number of the respective estimated segments 
K and the replicas' information theoretic measures (normalized mutual information 7 N and variation of 
information V) as a function of the resolution parameter 7. Decreasing 7, the minima of Eq. (1) lead 
to solutions progressively lower intra-community edge densities, effectively "zooming out" toward larger 
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structures. We determine all natural network resolutions by identifying the values of 7 for which the replicas 
exhibit extrema in the average information theoretic overlaps when expressed as a function of 7 [9]. For 
more detail sec [1]. 

3.2 Mixture of Gaussian Distribution Method 

The MGD based segmentation method models the data to be a mixture of multivariate Gaussian distributions 
of unknown means and covariances [5]. For 2D FLIM data, FLTs of all the pixels are assumed to be a 
mixture of univariate Gaussian distributions of unknown means and variances, and it is segmented using the 
expectation maximization (EM) based optimization method for segmenting mixture of Gaussian distributions 
developed originally by Hastie et al. [5] . We discuss below the MGD method for the general case of mixture 
of multivariate Gaussian distributions, assuming that the number of models is known. We then present how 
we propose using this method specifically for the 2D data as described in Section 2. 

The MGD method models the probability density function (pdf) of the fc th (V k € {1,2,... K}) component 
of the data x (x E M. d ) as Eq. (5), where K is the number of models defined by the number of Gaussian 
distributions in the mixture, d is the dimension of the data, and fi k of size d x 1 and of size dx d are the 
mean and covariance of the k th component. Assuming the weight of the k th component as a k , the mixture 
pdf is given by, Eq. (4). 



To estimate the segment identity for the observed incomplete data {x Xl x 2 , ■ ■ ■ , x N }, where N is the sam- 
ple size, we define complete data for performing the EM based optimization as {(xi,yi), (x 2 ,y2), ■ ■ ■ , {xn,Un)} 
Here yi (Vj/j e {1,2,..., K}) is the unknown segment identity of sample (Vi G {1,2,..., N}). The other 
unknown parameters 6 includes a*,, fi k , and ( k e {1,2,..., K}). The maximum likelihood function is 
given by Eq. (6), 



K 




(4) 



where 




(5) 




(6) 
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A classification EM algorithm developed by Hastie et al. maximizes Eq. (6). To segment a 2D NIR 
FLIM image, we segment the collection of the FLTs {x\ — T\,xi — r 2 , . . . ,xjy — t~n} from all the pixels 
of this image using the univariate version (d = 1) of the MGD algorithm. For simplicity, we assume K to 
be known, and perform the segmentation for increasing number of segments K. The final segmented image 
is generated using the K for which each of the resulting segments contain at least a predefined number of 
pixels, see Section 5. 

4 Segmentation of Images Described by 3D Spatio- Temporal In- 
formation 

This section extends the methods developed in Section 3 for spatially segmenting images described by 3D 
spatio-temporal information. 

4.1 Community Detection Method 

For spatially segmenting images described by 3D spatio-temporal signatures, we propose to apply multireso- 
lution CD method in two steps. In the first step, the multiresolution CD method identifies strongly-correlated 
replicas in the images based on their temporal signature, and finds hidden spatial structures in them. Here 
if a pixel pair is always in the same community (segment) in all the replicas, they must have a strong pref- 
erence to be part of the same segment or have a large edge weight. Similarly, if a pixel pair is not always 
in the same segment in all the replicas, they must have a preference not to be part of the same segment in 
all replicas or have a small edge weight. The resulting edge weights are in turn used in the second step to 
perform another multiresolution CD to spatially segment the target image. This strategy incorporates both 
the intensity and FLT information to determine spatial structures automatically in 3D FLIM data described 
by a spatio-temporal information. 

To differentiate the proposed approach with the general image segmentation method developed in our 
earlier work [1], we note that the method described in [1] first generates R replicas by permuting a 'symmetric' 
initial state defined by one pixel per segment of the studied system. It then applies the CD algorithm to 
each replica, and records the segment membership for each pixel. In contrast, for the images described by 
3D spatio-temporal signatures in this manuscript, the multiresolution CD method described in Section 3.1 
considers the spatial 2D frames of this 3D data as series of replicas, and segments each of them in its first 
step. Namely, each replica is represented by a frame corresponding to each temporal location of the 3D 
FLT data. The weight of each pixel pair is calculated based on the statistics of replicas, details of which are 
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described below. 

For performing the multiresolution CD in the second step, we define the edge weight Pij between the i th 
and j th € {1,2,..., N}) pixels of the input 3D image as follows, 



r=l 



where 



S a r a r = < 



1, if pixel i and j belong to the same segment in replica r, 

(8) 

0, otherwise. 



Using Eq. (7), we redefine the analogous Hamiltonian of Eq. (1) as, 
1 K \ 1 



which we minimize for performing the multiresolution CD. 
4.2 Mixture of Gaussian Distribution Method 

For performing the MGD based segmentation using the 3D spatio-temporal data, we normalize each temporal 
data corresponding to pixel i (Vi € {1,2,..., N}) in 2D spatial location to form Xj = x i2 , ■ . ■ , x i( i] T , 
where d = 256 in our case, see Section 2. The normalization here ensures max {xu\ I £ {1,2, ... = 1. 
The resulting Xj (i & {1,2, ... , N}) are used to perform the multivariate MGD as discussed in Section 3.2. 



5 Results 

This section describes representative examples of the proposed multiresolution CD method using an ex vivo 
liver tissue sample treated using an NIR fluorescent dye. We also compare the performance of our proposed 
method with the MGD method. 



5.1 Multiresolution Community Detection for Varying 7 

For the FLIM image of the liver tissue shown in Fig. 1A, we define the edge weight between two pixels as the 
absolute FLT difference between them. We applied the multiresolution CD to solve the resulting network 
formed by the image pixels as nodes. Fig. IB shows the plots of the information theoretic overlaps between 
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the replicas of the multircsolution CD, such as their normalized mutual information In and variation of 
information V, together with the number of estimated segments if as a function of the resolution parameter 
7. Decreasing 7, the minima of Eq. (1) leads to solutions with progressively lower intra-segment edge 
densities, effectively "zooming out" toward larger structures. We determine all natural network resolutions 
by identifying the values of 7 for which the replicas exhibit extrema in the average of their information 
theoretic overlaps when expressed as a function of 7. The extrema and plateau of the average information 
theoretic overlaps as a function of 7 over all replica pairs indicate the natural network resolutions [9]. 

Figs. 1C-1G show the image segmentation results by our multiresolution CD algorithm at different res- 
olutions. As the resolution decreases from Fig. 1C to 1G, the images show less detailed structures. In Fig. 
1C, three segments are clearly visible (light blue, green, and orange), in addition to the background segment 
(brown). In Fig. 1G, only one segment (light orange) appears against the background (light blue). Figs. 
1D-1F show the results produced using the multiresolution CD algorithm with the ranges of 7 that are in 
between the ones used for generating Figs. 1C and 1G. The segmented images show two distint segments 
located in the background. Thus by using different resolutions 7, the multiresolution CD method is able to 
detect the structures at different scales. 

5.2 Performance for Images Described by 2D Spatial Information 

We then compare in Fig. 2 the multircsolution CD method with the MGD method for the 2D spatial data 
of the FLIM image shown in Fig. 1A. In Figs. 2A and 2C, we present the image segmentation results 
by the MGD and multiresolution CD, respectively. Note that we evaluated the former method here for 
increasing number of segments K, and for K > 3, the method started producing segments with pixels fewer 
than 100. In contrast, the multiresolution CD method selects the number of segments automatically as 
discussed in Sections 3.1 and 4.1. In the current example, we used 7 € (0.6, 0.7) for performing segmentation 
using this method. To measure the performances of the two methods, we plot in Figs. 2B and 2D the 
respective normalized average decay traces as the function of photon arrival time in the PMT. The proposed 
multiresolution CD method is able to identify segments with more distinct and diverse decay traces than the 
other method as seen in these two plots. This method allows the user to adjust the resolution 7, and thus 
to obtain the segments with more distinct and diverse decay traces. 

To quantify the diversity offered by the respective segmentation methods, we compute the Shannon index 
for each of them. Shannon index has been widely used in the literature to quantify species diversity in a 
habitat or community [16]. We compute this index for our study using the set described by the Euclidean 
distances between every pairs of average decay traces of the estimated segments for MGD or multircsolution 
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CD method. To perform this computation, we first bin the computed distances using intervals defined by 
(0.2(m - 1), 0.2m], where m € {1,2, ... , M} and 



M = 



max 



(c^mgd, c^mcd) 



(10) 



0.2 



where ["•] denotes a ceiling function and c^mgd and c?mcd are the maximum computed Euclidean distances 
between the average decay traces of the estimated segments for MGD and multircsolution CD methods, 
respectively. We then compute the frequency of occurance of the distances in each bin, and denote this 
frequency by q m for the m th bin. We finally compute the Shannon index for MGD or multiresolution CD 
method as, 



We obtained Hmgd = 0.64 and -S/mcd = 1-26 for the MGD based and multiresolution CD based segment 
estimates of the 2D spatial data of the FLIM image shown in Fig. 1A. This result suggests that the mul- 
tiresolution CD method offers almost more than two times diversity than the MGD method in segmenting 
FLIM images described by 2D spatial information. 

5.3 Performance for Images Described by 3D Spatio-Temporal Information 

We finally compare in Fig. 3 the multiresolution CD with the MGD method for the 3D spatio-temporal 
version of the FLIM image shown in Fig. 1A. In Figs. 3A and 3C, we present the image segmentation 
results by the MGD and multiresolution CD, respectively. The optimal number of segments was selected 
respectively by the MGD method and the CD method using the similar procedures as stated in Sections 3 
and 4. To conduct the segmentation using the multiresolution CD on the 3D spatio-temporal data, we used 
single resolutions in both steps; particularly, we used 71 = 1 and V\ = 2.5 in its first step and 72 = 10 and 
V2 = 2.5 in its second step. Both of the methods are able to produce more connected structures for the 3D 
spatio-temporal data than was obtained for the 2D spatial data. The resulting normalized average decay 
traces of the estimated segments are shown in Fig. 3B and 3D for the respective methods. For a clearer 
depiction, Fig. 3D here shows the decay traces only for segments with number of pixels more than 300. 
To quantify the diversity offered by the respective segmentation methods, we computed Hmgd = 1.1 and 
Hmgd = 1-53 for the MGD based and multiresolution CD based segment estimates of the 3D spatio-temporal 
version of the FLIM image shown in Fig. 1A. This result suggests that the multiresolution CD method offers 
almost more than 1.25 times diversity than the MGD method in segmenting FLIM images described by 3D 



M 



H = - ^2 9m log 9m- 



(11) 



771—1 
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spatio-temporal information. 

6 Conclusion 

We have developed a multiresolution community detection algorithm to segment two-dimensional and three- 
dimensional fluorescence lifetime imaging microscopy (FLIM) data. The proposed method is able to identify 
structures in different scales in the input FLIM images based on an information theoretic measures. It offers 
atleast 1.25 times more diversity based on Shannon index than a state-of-the art method based on mixture 
of Gaussian distributions, in producing distinct and diverse decay traces in the segmented images. 
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A Definitions: Trials and Replicas 

We review the notions of "trials" and "replicas," used in our community detection (CD) algorithms. Both of 
these notions pertain to the use of multiple identical copies of the same system which differ from one another 
by a permutation of the initial site indices. Thus, whenever the time evolution depends on sequentially 
ordered searches for energy lowering moves (as it does in our greedy algorithm), these copies may generally 
reach different local solutions. By the use of an ensemble of such identical copies, we attain accurate results 
as well as determine information theoretic correlations between the candidate solutions, and infer from them 
a detailed picture of the system. 

In the definitions of "trials" and "replicas" given below, any given algorithm may be used to minimize 
the selected cost function. In our particular case, we minimize the Hamiltonian of Eq. 1. 

• Trials: We use "trials" alone in our bare community detection algorithm. We run the algorithm on the 
same problem T independent times. This may generally lead to different contending states that minimize 
Eq. 1. Out of these T trials, we pick the lowest energy state and use that state as the solution. 

• Replicas: We use both "trials" and "replicas" in our multiresolution CD algorithm. Each sequence of 
the above described T trials is termed as a replica. When using "replicas" in the current context, we run 
the aforementioned T trials (and pick the solution that attains lowest energy in the Hamiltonian of Eq. 



12 



1) R independent times. By examining information theoretic correlations between the R replicas, we infer 
which features of the contending solutions are well agreed on (and thus are likely to be correct), and on 
which features there is a large variance between the disparate contending solutions that may generally mark 
important physical boundaries. We compute the information theoretic correlations within the ensemble of R 
replicas. Specifically, the information theoretic extrema as a function of the resolution parameter, generally 
correspond to more pertinent solutions that are locally stable to a continuous change of scale. It is in this 
way we detect the important physical scales in the system. 

B Information Theoretic Measures 

We use information theoretic measures to calculate correlations between community detection (CD) solutions. 
The CD method partitions N pixels for a replica r (Vr G {1, 2, . . . , R}) into K r segments, where segment k 
(k G {1, 2, . . . , K r }) consists of pixels. The ratio N^/N is the probability that a randomly selected pixel 
is found in the segment k (k G {1,2,..., K r }). The Shannon entropy [11] is 



where Nk 1 k 2 is the number of common pixels in the segment k\ {k\ G {1, 2, . . . , K r }) of replica r (r G 
{1,2,..., R}) and the segment k 2 (k 2 G {1, 2, . . . , K s }) of replica s (sg {1,2,..., R}). 
The variation of information V(r, s) between the two segments r and s is 




(12) 



The mutual information I(r, s) between the replicas r and s ({r, s} & {1,2, ... , R}) is 




(13) 



V{r,s)=H r + H s -2I{r,s), 



(14) 



which has a range of < V(r, s) < log 2 N. 

The normalized mutual information I^s(r, s) is 



Jn(t",s) = 



2J(r, s) 
H r + H s ' 



(15) 



with the obvious range of < I^(r, s) < 1. 



13 



Higher In(-) and lower V(-) values indicate better agreement between the compared segments. 

C Community Detection Algorithm 

Our community detection (CD) algorithm for minimizing Eq. (1) follows four steps [10]. 

1. We partition the pixels based on a symmetric or fixed K initialization where K is the number of 
segments. 

• Symmetric initialization alludes to an initialization, where each pixel forms its own segment; i.e., 
initially, there are = N segments. 

• Fixed K initialization corresponds to a random initial distribution of all pixels into K segments. 

For image segmentation, a symmetric initialization is used for the unsupervised case. In this case, the 
algorithm docs not know what the number of segments are, so the symmetric initialization provides 
the advantage of no bias towards a particular segment. The algorithm decides the number of segments 
K, by means of the lowest energy solution. 

Fixed K initialization may be used in a supervised image segmentation. The community membership 
of an individual pixel is then changed to lower the solution energy using CD algorithm. Here the user 
has to decide how much information is needed by observing the original image, and enter the number 
of segments K as an input, where different levels of information correspond to different numbers of 
segments K. For instance, if only one target needs to be identified, K = 2 is enough, which describe 
the target and background. 

In our work, we use the unsupervised image segmentation, and let the algorithm decide the community 
number K. 

2. Next, we select each pixel, and place it in the segment that best lowers the energy of Eq. (1) based on 
the current state of the system. 

3. We repeat this process for all pixels, and continue iterating until no energy lowering moves are found 
after one full cycle through all pixels. 

4. We repeat these processes for T trials, and select the lowest energy result as the best solution. Different 
trials differ solely by the permuted pixel order of the initial state. In the multi-resolution application 
discussed below, we further use different replicas, i.e., different independent solutions attained using 
our proposed algorithm, and solve the system based on a finite temperature algorithm [11]. 
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D Mult iresolut ion Algorithm 

We illustrate below how the multiresolution algorithm [9] works. 

To begin the multiresolution algorithm, we need to specify the number of replicas R at each resolution, 
the number of trials per replica T, and the starting and ending resolutions, 70 and 7/, respectively. We 
typically use the number of replicas as 8 < R < 12, and the number of trials as 2 < T < 20. We select the 
lowest energy solution among the T trials for each replica. The initial state of the replicas are generated by 
permuting the pixel labels in the symmetric initialized state of one pixel per community. These permutations 
P simply reorder the pixel indices (1, 2, 3, . . . , i, . . . , N) — > (PI, P2, . . . , PN) (with Pi the state of i under a 
permutation), and thus lead to a different initial state. 

1. The algorithm starts from the initialization of the system, described in item (1) of Section C. 

2. We then minimize Eq. (1) independently for all replicas at a resolution 7 = 7$ € {70, 71, ... , 7/-i,7/} 
as described in Section C. Initially, i = (i.e., 7 = 70). 

3. The algorithm then calculates the average inter-replica information theoretic measures, such as Jn and 
V, at that value of 7. 

4. The algorithm then proceeds to the next resolution -f i+ i e {70,71, ■ • • ,7/-i,7/} (with -f i+ i > 7,). 

5. We then return to Step 2. 

6. After examining the case of 7 — jf, the algorithm outputs the inter-replica information theory overlaps 
for the entire range of the resolutions studied. 

7. We examine those values of 7 corresponding to the extrema in the average inter-replica information 
theoratic overlaps. Physically, the resulting image segmentation for these values is locally insensitive 
to the change of resolution (i.e., small changes in 7) and generally highlights prominent features of the 
image. 

With A and B denoting graph segments in two different replicas, and Q(A, B) their overlap, the average 
inter-replica overlap for a general quantity Q(-) [9] is explicitly 




(16) 



Similarly, for a single replica quantity, the average is, trivially, (Q) = J2aH(^)/R [9]. 
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E Classification Expectation Maximization Algorithm 

To maximize Eq. (6), we use a classification expectation maximization algorithm [5], steps of which we briefly 
review below. For more discussion, readers might consult the classical book on machine learning written by 
Hastie et al. [5]. 

1. We initialize the parameters a k , and (k G {1, 2, . . . , K}). To perform this initialization, we 
first initialize the segment identity tji (Vt/j £ {1,2, ... , K}, i e {1,2,..., N}) by randomly assigning its 
value to be 1, 2, . . . , or K. We perform this random assigment here by uniformly drawing a number 
from {1,2,..., K}, and assign this number to yf * 1 (i e {1,2,..., N}). Based on this assignment, we 
compute a k °^ (k € {1, 2, ... , K}) as, 

jo) _ Eti 1 (vf ] = k ) 

a k — ~jy i \ L <) 

where 1 (y^ = k) is an indicator function as given by, 

nyf } = k) = l (18) 

0, otherwise. 
Similarly, we compute /zj^ and (fee {1,2,..., K}) as, 

„(o) _ EL 1 (y? ] = k ) x i nq x 
Mfc " EL^ 0) = *) ' ( 9) 



and 

2. Expectation step: This step computes the posterior probabilities for all {i = 1,2,..., N} and {k = 
1,2,..., K} for the (/ + l) th iteration, using updated estimates of the parameters a k , and from 
the I th iteration. This step determines how likely Xi is to be the part of the k th segment. The posterior 
probability is given by, 



a^Ux \u {l) S (() ) 

_ _ a k Jk\Xi\V k '^fc ) / 91 \ 
Pl - k ~ (I) , , | (!) ^(Ox' 1 J 



1G 



3. Classification step: This step determines for the (I + l) th step which among k = 1, 2, . . . , K is most 
likely to be the segment identity of Xi (i & {1,2, ... , N}). Namely, we compute, 



J> = argmax p iik , 



(22) 



Pi,k> 



1, if fe' = argmax Pi t k, 
k 

0, otherwise. 



(23) 



4. Maximization step: Using the information obtained from Eqs. (16) and (17), we update the estimates 
of the parameters Ofc, /x fe , and as, 



(1+1) _ L,j=l Pi,k 

TV ' 



(24) 



v-^AT „ 
„(«+l) _ 2^i=lPi,kXi 
— - 



(25) 



(26) 



5. We repeat Steps 2-4, until convergence. The convergence is determined by 



|L({-}|e (,) ,»«)-L({-}|fl (, - 1) ,y( , - 1 ))| 



< c, 



(27) 



where e is an infinitcsimally small number, which was 0.001 in our case. The EM algorithm described 
herein monotonically increases the cost function defined by Eq. (6) [5], and the unknown parameters 
are guaranteed to converge, at least locally. 
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Figure 1: Segmentation results for an ex vivo liver tissue image with two-dimensional (2D) spatial infor- 
mation: Panel A shows the result generated by the multiresolution community detection (CD) method for 
segmenting images described by the 2D fluorescence lifetime imaging microscopy data of an ex vivo liver 
tissue sample. Panel B shows the normalized mutual information In and the variation of information V of 
the replicas of the multiresolution CD method and its number of estimated segments K as the function of 
resolution parameter 7. We selected the values of 7 at the peaks of the curve described using V, and obtained 
the segmented images as shown in Panels C-G by our proposed multiresolution CD method. Images with 
the higher values of 7 have more precise structures, see Panels C and D. Images with the lower values of 7 
have rough structures, see Panels E, F, and G. 
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Figure 2: Comparison between the multiresolution community detection (CD) method and the mixture of 
Gaussian distributions (MGD) method for images described by two-dimensional spatial information: Panel 
A is the segmentation result by the MGD method. Panel B shows the normalized average decay traces 
corresponding to the estimated segments shown in Panel A. Panel C is the segmentation result by the 
multiresolution CD method using 7 = 0.7. Panel D shows the normalized average decay traces corresponding 
to the estimated segments shown in Panel C. Multiresolution CD is able to identify more distinct and diverse 
decay traces than the MGD method. 



21 





Photon arrival time (ns) 




Figure 3: Comparison between the multiresolution community detection (CD) method inspired by the replica 
correlation and the mixture of Gaussian distributions (MGD) method for images described by three dimen- 
sional (3D) spatio-temporal information. Panel A is the result obtained by the MGD method. Panel B shows 
the normalized average decay traces corresponding to the estimated segments shown in Panel A. Panel C 
is the segmentation result by the multiresolution CD method developed for images described by 3D spatio- 
temporal information, using 71 = 1 and V\ = 2.5 in its first step and 72 = 10 and V2 — 2.5 in its second 
step. Panel D shows the normalized average decay traces corresponding to the estimated segments shown 
in Panel C. Both of the methods are able to produce more connected structures for the 3D spatio-temporal 
data than that was obtained for the 2D spatial data. Our results indicate that both of the methods here 
appear to be performing similarly in terms of generating segments with distinct and diverse decay traces. 
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